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ABSTRACT 



An adaptation of a primitive variable, finite-difference 
computer program was accomplished in order to predict the 
non-reacting flow fields in turbojet test cells and the 
reacting flow fields in solid fuel ramjets. The study 
compares the predictions of the primitive variable computer 

model with an earlier computer model and empirical data. 

\ 

It was found that the new model reasonably predicted the 
flow fields in both geometries. In addition, the primitive 
variable model allowed simulation of test cell flows up to 
full engine throttle conditions and solid fuel ramjet flows 
which included an aft mixing chamber. 
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ABBREVIATIONS AND SYMBOLS 



ROMAN 
A 
BP 
C- 



SYMBOLS 

Area 

Mass transfer (or "blowing") parameter 
K-e model empirical constants (Table I) 



E 

g 

G 

gr 

h 

h 

H 

h 

i 

I 

K, k 
K 
1 
m 
M 
M 
m 
m " 



Coefficient of friction 

Specific heat at constant pressure 

9.0 

Mass transfer conductance 

Air mass flux 

Gram 

Enthalpy 

Stagnation enthalpy 
Dimensionless enthalpy 
Heat transfer conductance 
Stoichiometric coefficient 
Turbulence intensity 
Turbulence kinetic energy 
Thermal conductivity 
Length scale of turbulence 
Mass fraction 
Molecular weight 
Average molecular weight 
Mass flow 
Mass flux 
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Pressure 



P,P 

q" Heat flux 

r Radial distance 

R Gas constant 

R Universal gas constant 

RR, r Regression rate 
S Source terms 

St Stanton number 

T Temperature 

u Axial velocity 

v Radial velocity 

x Axial distance 

y^ + Dimensionless distance from solid boundary 

GREEK SYMBOLS 



r 

5 

AH 

£ 

K 

u 

p 

a 

T 

X 



Effective transport coefficient 

Incremental distance from wall 

Heat of combustion per Kg of fuel 

Turbulence dissipation rate 

von Karman constant 

Viscosity 

Density 

Prandtl or Schmidt Number 

Shear stress 

Any variable 

itu - m /i 
fu ox 7 
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air Air 

atm Atmosphere 

bw Fuel surface (or "blowing wall") 

c Conserved 

eff Effective 

fg Fuel grain 

fp Fuel port 

fu Fuel 

in Inlet 

lam Laminar 

N2 Nitrogen 
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p Near wall node 
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ref Reference 
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I . INTRODUCTION 



A. BACKGROUND 

During the past few years, there have been many advance- 
ments in numerical techniques for predicting the behavior 
of fluid flows. For example, several computer models have 
been developed by Gosman, Spalding and others [1,2,3] which 
use the mass, momentum and energy conservation equations 
reduced to finite difference nonlinear algebraic form. The 
development of reliable computer programs of this type 
greatly benefits engineering analysis in such widely varying 
fields as meteorology, aerodynamics and gasdynamics. 

The earlier two-dimensional computer codes were based on 
vorticity (oj) and stream function (ip) [1,2,5]. This form 
of the governing equations eliminates pressure and velocity 
from immediate consideration. Pressure is normally calcu- 
lated only after a converged solution is obtained. This 
technique has several inherent disadvantages: 

1. It results in large errors in the predicted pressure 
distributions in all but quiescent flow regions due to the 
higher order dependence of the pressure gradient on stream 
function [6]. 

2. It is usually restricted to constant density flows 
or to flows in which density varies only with temperature 
[3,6] . 

3. The boundary conditions are difficult to specify 
[3,4,5] . 
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4. Considerable difficulty has been experienced in 
arriving at converged solutions, especially for non- 
uniformly spaced grids and high flow rates [2, 5, 6, 7]. 

5. The \p-Ci) model is not easily extended to three 
dimensional flows [3] . 

To overcome these difficulties, emphasis has been placed on 
developing computer codes based on velocity and pressure, 
the primitive variables. 

A major problem with any new computer model is model 
validation. The difficulties of collecting accurate empiri- 
cal data are multiplied when investigating three dimensional 
and/or reacting flows. In addition, many variables within 
these flows are not readily measurable (turbulence intensi- 
ties , etc . ) . 

An effort to utilize elliptic computer models which 
can handle turbulent, reacting, variable density flows at 
high subsonic and sonic velocities has been underway at the 
Naval Postgraduate School for several years. Two specific 
areas which have been investigated are flows in a turbojet 
test cell and in the combustion environment of a solid fuel 
ramjet. 

B. TURBOJET TEST CELL 

It is important to have the capability to test high 
performance jet engines throughout their operating envelope 
under conditions which approximate installed conditions. 

This is accomplished in blockhouse type installation called 
turbojet test cells (TJTC) . The typical test cell incorporates 
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an inlet, a horizontal test section and a vertical exhaust 
stack. The engine to be tested is normally mounted near 
the center of the cell to allow the development of a nearly 
uniform engine inlet velocity profile. The engine exhausts 
into an augmentor tube which entrains additional air for 
exhaust gas cooling and dilution. The quantity of this 
secondary air is crucial to proper engine testing and test 
cell performance. Each cell is equipped with adequate instru- 
mentation to allow assessment of engine performance parameters. 

Testing today's high power and high mass flow engines in 
these installations produces a myriad of noise and air pollu- 
tion problems. Cell modifications must often be made to 
minimize these problems. This fact coupled with the future 
need for larger, more expensive test cells to replace obso- 
lete cells and to accommodate new generations of high tech- 
nology engines, makes the development of reliable modeling 
methods imperative. The frequently used one-dimensional models 
are not adequate for predicting the details of the compli- 
cated flows within a turbojet test cell and, therefore, 
the cells often do not perform to their designed limits. 

An accurate flow model would provide a needed design tool 
which could help prevent costly design errors and improve 
operating efficiency. 

A two-dimensional tp-u computer code was used to analyze 
the flows in a full scale and a subscale turbojet test cell 
at the Naval Postgraduate School [2,6]. Experimental data 
from the subscale test cell have been compared with 
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computations made with the computer model [ 6 J and the latter 
has aided in design modifications and the evaluation of 
pollution control equipment on that installation. However, 
as discussed above, this analysis technique has several 
disadvantages . 

A primitive variable (u-v-p) computer model could in- 
crease this prediction capability by extending it to specific 
geometries and flow rates that the \p-u model is incapable of 
predicting. In addition, a u-v-p model would more readily 
allow variable density flows to be analyzed and should more 
accurately predict augmentor pressure distribution. 

C. SOLID FUEL RAMJET 

A solid fuel ramjet (SFRJ) most often consists of a solid 
fuel grain which provides the walls for the combustion cham- 
ber [5] . Located at the air inlet end of the combustor is 
a sudden expansion or other type of flame stabilization 
device. The opposite end, downstream of the fuel grain, 
may also incorporate a sudden expansion aft mixing chamber. 
The primary combustion region is a turbulent diffusion flame 
which emanates from the forward recirculation zone and 
remains within the developing boundary layer. The aft 
mixing region may incorporate some means of injecting air 
for burning the fuel-rich mixture which has been found to 
exist there [5], Mixing chamber and inlet design variables, 
fuel grain design and fuel properties make a wide variety 
of performance characteristics available. 
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The possibility of incorporating this type of propulsion 
device into future medium or long-range tactical weapon 
systems coupled with the expense of testing each new design, 
makes the development of a reliable computer model of this 
system highly desirable. The model could be used to pre- 
dict the effects of fuel properties and to inexpensively 
evaluate different geometries and operating conditions. In 
addition, a three dimensional code would allow modeling dis- 
crete air injection into the aft mixing region. The latter 
technique can substantially increase combustion efficiency 
and allowable fuel loading. 

Previous work at the Naval Postgraduate School has been 
directed toward improvement of the quantitative accuracy 
of the ip— to model and toward validation of that model [5] . 
Reasonable agreement with empirical data has been obtained. 
However, as previously stated, the ip-oj model does not pre- 
dict accurate pressure distributions and numerical diffi- 
culties prevented modeling the aft mixing chamber. 

The purpose of this investigation was to adapt and 
validate a primitive variable, two-dimensional, finite 
difference computer code which models the flows within 
turbojet test cells and solid fuel ramjets. 
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II. MODEL OVERVIEW 



A. INTRODUCTION 

The computer model used in this study was adapted from 
the CHAMPION 2/E/FIX computer program developed by Pun and 
Spalding. As described in detail in reference 8, CHAMPION 
is a TWO- dimens ional Elliptic, Fix ed grid computer program 
which provides a solution of the conservation equations for 
recirculating flows in finite difference form. 

B. ASSUMPTIONS 

The flow was assumed to be steady, two-dimensional and 
subsonic. For simplicity the value of specific heat (C ) 

hr 

was assumed to be constant although its dependence on 
temperature and/or composition could easily be included. 

A modified Jones-Launder [8,10,11,12] two parameter 
turbulence model was incorporated to calculate the effec- 
tive viscosity. It uses five empirical constants (Table I) 
and requires that two additional variables , turbulent kinetic 
energy (K) and turbulent dissipation rate (e) , be evaluated. 
Effective viscosity was calculated using the formulas: 

p ef f = p lam + “t (1) 

where 

P t = C D p K 2 /e (2) 
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For reacting flows, the four species, oxygen, nitrogen, 
fuel, and products, were considered. Simple, one-step, 
infinitely fast kinetics were assumed in which a fuel com- 
bines with an oxidant to form a single product without 
intermediaries [5,9]. 



1 gr fuel + i gr oxidizer -> (1 + i) gr products 



Fuel and oxygen, therefore, could not exist simultaneously 
and the combustion process was mixing limited. In addition, 
it was assumed that no oxygen existed at the fuel surface 
and that surface was isothermal. The turbulent Prandtl and 
Schmidt numbers were taken equal to unity and, therefore, 
the turbulent Lewis number was unity. The laminar Prandtl 
number was also taken to be unity. 

C. GOVERNING EQUATIONS 

The conservation equations for axi- symmetrical flows 
with no tangential variations can be put into the general 
form [8] : 



3 , . 13 , , 

33j( p u*) + ? jjlprv*) 

convection terms 



- i_ (r 3<K _ 1 i_ 
3x U d)3x 



Sr' r % 



it) 

3r 



diffusion terms 



S A 

<p 

source 

terms 



where <j> stands for the dependent variable (u,v,k,e,h, etc.) 

being considered (j> = 1 for the continuity equation) , r is 

9 

the appropriate effective exchange coefficient for turbulent 
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flow and S A is the "source term" (Table II) . The energy 
equation in terms of stagnation enthalpy has no source terms 
since the turbulent Prandtl and Schmidt numbers were chosen 
as unity and radiative transport was neglected [1,3]. Thus 
the stagnation enthalpy is given by: 

h = h + (u 2 + v 2 )/2 + K (4) 

where for non-reacting flows: 

h = C T (5) 

P 



and for reacting flows: 

h = rn AH/i + C (T - T -) (6) 

ox p ref 

The calculation of temperature was made using equations 
(4) , (5) and (6) . Density was calculated from the perfect 

gas law: 



P/RT 




(7) 


for non-reacting flows: 


R = constant 


(8) 


for reacting flows: R = 


R/M 





where , 



l/M = "i fg /M fg + i» ox /« ox + ”n 2 /M N2 + m pr /M pr 

(9) 
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For modeling reacting flows, two additional quantities, 
and X = - m /i, were evaluated. Each of these 

properties as well as stagnation enthalpy have identical 
governing differential equations (equation (3) with no 
source terms) . In appropriate dimensionless form they also 
have identical boundary conditions. Thus, only one of the 
equations had to be solved. The dimensionless form selected 
for each property was: 



H 



< h in - h)/(h in - V 



( 10 ) 



m N2 (m N2- m N2 ) / (m N2.“ m N2. } 

in in rg 

X = (X - X in )/(x £g - x ln > 



In this study, stagnation enthalpy was calculated. H was 
then formed using equation (10) . Since H = m^ = x" at all 
points in the flow field, m^ and X could be calculated 
using equations (11) and (12). The mass fractions of fuel, 
oxygen, and products ( m £ u / m ox ' m pr^ were found from the 
equations : 



for x > 0; m- = X/ m 

A — f u A ox 

for x < 0 ; m- = 0 , m 

A fu ox 



0 



(13) 

“X i 



m 



Pt 



= 1 - m 



ox 



- m 



fu 



^2 



(14) 
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D. CONSERVATION OF MASS 



On each radial line the mass flow rate was calculated 
using the local density. The error in mass flow (compared 
to the summation of "mass-in" at all upstream boundaries) 
was used to uniformly adjust the axial velocity over the 
entire line. This process ensured that overall continuity 
was satisfied on the line. The pressure at all downstream 
locations was then adjusted to approximately correct for the 
momentum imbalance created by the uniform axial velocity 
adjustment. A "pressure correction" equation was then 
solved for each cell on the line. Local cell velocity 
(axial and radial) and pressure were then adjusted to satisfy 
cell-wise continuity. The details of this procedure are 
presented in reference 8 . 

E. BOUNDARY CONDITIONS 
1 . Introduction 

Fixed boundary conditions were set at the desired 
or experimentally determined value and held constant. 
Specified gradient boundary conditions were handled by 
setting the appropriate convection/diffusion coefficient 
to zero in the finite difference equation ("breaking the 
link") and then entering the appropriate gradient through 
linearized "false" source terms [8] . The geometries and 
appropriate boundary conditions for the TJTC and SFRJ are 
summarized in figures l and 2. 
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2. Inlet 



Although not a computer program limitation, "plug 



flow" was assumed at the inlets for both the SFRJ and TJTC. 
The secondary flow inlet of the TJTC was recessed approxi- 
mately 0.3 meters (figure 1) to allow a velocity profile to 
develop over the length of the engine. Turbulent kinetic 
energy was selected to be uniform with a value which corres- 
ponded to the approximate turbulent intensity of the inlet 
flow. 

3 . Axis of Symmetry and Exit Plane 



on the center line and exit respectively, with one exception 
the radial flow velocity was zero. 

4. Solid Boundaries 



adiabatic with both velocity components equal to zero ("no 
slip" condition) . 



The border between the laminar sublayer and the turbulent 



Radial and axial gradients were set equal to zero 



All non-reacting solid boundaries were considered 



For simplicity, a two part boundary layer was used. 



layer was taken at y + = 11.5 [8] . y ' was evaluated at 

P P 

each near wall node (p) , 



y. 



P 




(15) 



where, for y. 



+ 



> 11.5 



P 



T 



W 




(16) 
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was assumed uniform from the wall to the near wall 



w 



grid point. Thus, 



+ _ 1/4 __ 1/2 . . 

= c d p K P 6/u 



lam 



( 17 ) 



If y > 11.5, the wall shear stress (x ) was calculated 
u p — w 

using the formula: 



T w - C D 1/2 p K p = f C D 1/4 K p 1/2 (u/u + ) 



1/4 



1/2 



- k C„"’ p u p V' Vln<EoSC D 1 / 4 K p 1 / 2/,J lain ) 



(18) 



where 



u + = — ln(E y + ) 

< J p 



(19) 



Wall shear stress was evaluated for 
formula: 



+ 



< 



11.5 from the 



: = u n u /6 

w lam p 



( 20 ) 



Due to the steep gradients of properties in turbulent 
flows near solid boundaries, the source terms for K and e at 
near wall nodes were expressed in terms of the wall shear 
stress [1,8]. x w also provides the boundary condition for 
the u and v equations. In the following equation for 
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turbulence dissipation rate (e) at a near wall node (p) , 
the length scale is presumed proportional to the distance 
from the wall (6) . 

_ 3/4 3/2. ~ „ 3/2 

£ p ~ C D K p /<5 = K p /2 . 43 6 \21) 

It was found, as was previously found by Netzer [5] , that 

when using the sudden expansion geometry in reacting flows 

the near wall dissipation had to be increased on the step 
3 /2 

face (e = K / /0.4<5) and that the grid spacing adjacent 
P P 

to the fuel surface had to be fine (y + < 11.5) in order 

P 

to obtain a temperature distribution in qualitative agree- 
ment with experiment. Equation (20) implies that the wall 
shear stress is calculated assuming a linear velocity pro- 
file when y + < 11.5. A near-wall grid point, therefore, 
can lie within the laminar sublayer, but the source terms 
for K and £ imply that muc h greater than one 

[10,11]. This fact precludes y + from being significantly 

P 

less than 11.5. 

For reacting flows, the boundary conditions for the 
dimensionless properties (equations (10) , (11) and (12) ) 

were zero at the inlet and unity "deep" in the fuel grain 
(fg) . These properties were considered to have zero gradients 
on non-reacting surfaces. 

The assumptions employed for reacting flows (unity 
turbulent Prandtl and Schmidt numbers , simple chemical 
reaction, constant specific heat and stagnation enthalpy 
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defined in equations (4) and ( 6 ) ) result in a general 
boundary condition for all "conserved" properties ( 4 ^) 
[ 9 ] on a surface which has mass transfer, 



m 



It 

bw 



(r , inr’bv/^c 



bw 




( 22 ) 



where <j> c represents h, or X = m f u “ m /i. 

A mass transfer conductance (g) is often defined 
such that. 



(r 



dr 



bw 



= g U, 



'bw 



(23) 



where 4> is defined as the free steam value. For this 
c 

00 

application, <j> c was taken to be the local near wall value 




Substituting equation (23) into equation (22) 

yields : 




g(4>~ ~ ^ )/(<}> “ 4> ) (24) 

p bw c bw c f g 



= g BP 



(25) 



where BP represents the mass transfer (or "blowing") 
parameter. 

Without mass transfer the wall heat flux (q ") can 

w 

be defined in terms of the conditions at the near wall node. 
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It 



(26) 



q 



w 





, K ah ' 

X 3r 
P 



w 



where h w is the enthalpy of the wall and h is the heat 
transfer conductance. 

With 4> c = h in equation (23) , 

g = (5- 4^) /(h - h ) 

^ C 3r p w 

p w * 

Substituting equation (26) into equation (27) 

yields : 



g = 



ft/C 



P 



(28) 



or 



g/(pu) = Ji/[(pu) c] = St (29) 

thus , 

g = (pu) St (30) 

From Reynold's Analogy with unity Prandtl number, 
St = C f /2 = x w /(pu 2 ) p (31) 

where is the local friction coefficient. 
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Combining equations (30) and (31) yields: 



g 



t /u 
w p 



( 32 ) 



Using the Couette flow approximation for the 
boundary layer behavior with mass transfer [9] , 



g = g* In (1 + BP)/BP 



(33) 



where 



g* = lim Bp^ 0 (g) < 34 > 

In this application, BP was evaluated from the 
solution of the energy equation using, 

BP = (h - h, ) / (h, - h c ) (35) 

p bw /x bw fg 

The wall shear stress was calculated using equation 
(18) or equation (20) and modified with equation (33) . 

t bw = In (1 + BP) /BP (36) 

where x is the wall shear stress without wall mass addition, 
w 

The mass transfer conductance (g) was found using 
equation (32). The wall mass flux was then evaluated using 
equation (25) . 
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The wall heat flux (<3 W ") on all solid isothermal 
boundaries was evaluated using Reynold's Analogy. 

-q "/(h - h ) = x /u (37 

^w ' ' p w p 

q w " provides the boundary condition for the h equation. 

Since the blowing rates were small for the solid 
fuel ramjet (typically, BP < 2.0), K and e were evaluated 
using equation (3) and the terms presented in Table II 
which incorporate the empirical constants of Table I. 

Blowing velocity (v^) and fuel regression rate 
(RR) were calculated using the formulas: 



v, = -m, "/p, 

bw bw bw 



(38) 



RR = * bw "/p fg (39) 

F . SOLUTION PROCEDURE 

Five variables (u,v,K,e and H or h) were solved using 
equation (3) in finite difference form. The line by line 
iterative procedure employed upwind differencing and under 
relaxation to promote convergence [8] . Pressure (relative 
to a selectable position and magnitude within the grid) was 
obtained from the mass conservation imposed on each radial 
grid line and on each nodal control volume as discussed 
above. Effective viscosity, temperature and density were 
also obtained as described above. A more detailed explana- 
tion of this procedure can be found in reference (8). 
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III. RESULTS AND CONCLUSIONS 



A. TURBOJET TEST CELL 
1 . Introduction 

The purpose of this portion of the study was to 
utilize a primitive variable, finite difference computer 
program to analyze the flow within a turbojet test cell and 
to validate the model with data obtained from a subscale 
turbojet test cell located at the Naval Postgraduate School. 
Previous work [2,6,7] had accomplished this task for a 
computer model and, therefore, empirical data and the pre- 
dictions of the ip-co model were available for comparison. 

The experimental data available consisted of augmentor wall 
pressure distributions and radial velocity profiles along 
the length of the augmentor tube for low, medium, and high 
engine flow rates. 

As previously indicated, the model did a poor 

job of predicting pressure distributions in all but quiescent 
flow regions. Additionally, numerical convergence was diffi- 
cult to obtain with that model when used to predict high 
velocity flows where compressibility effects are significant. 
It was anticipated that a primitive variable model would 
help to alleviate these difficulties. It is desirable to 
have a model which can be used to predict the flow field for 
full-throttle engine conditions where the engine exhaust 
flow is choked. 
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Adaptation of the primitive variable model to the 
subscale test cell geometry required the use of several 
approximations : 

a. In modeling axi-symmetric flow, the engine was 
by necessity positioned at the axis of symmetry. In the 
actual test cell the engine was mounted closer to the deck 
than to the overhead. It would be expected, therefore, 
that the velocity distribution in the secondary flow (the 
flow around the engine) would be somewhat different than 
predicted. 

b. The subscale test cell cross section is rec- 
tangular while the engine and augmentor tube are cylindrical. 
The system was modeled as three concentric cylinders with 
cross-sectional areas equivalent to the physical system. The 
nozzle exit area, the test section cross-sectional area and 
the empirical augmentation ratios and mass flow rates were 
used in calculating the axial inlet velocities used in both 
models . 

c. The actual engine incorporated a converging noz- 
zle. The engine was modeled as a cylinder with a diameter 
equal to the actual nozzle exhaust diameter. 

d. In the model the augmentor inlet and the aft 
test cell wall were taken to be flush. The actual augmentor 
is often inserted into the test section which forms a recircu- 
lation zone above the augmentor tube. The effects on the 
augmentor flow field introduced by this recirculation region 
were shown to be minimal by Speakman, et al [7] . It should 
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be noted that the ip-u model incorporated this recirculation 
zone and a flow reducing lip flange on the augmentor inlet. 

When comparisons were made between the predictions of the 
two models, the effects of this recirculation zone and the 
augmentor inlet lip flange in the i|>-u> model were minimized 
by reducing their dimensions to one grid spacing. 

Three of the flow conditions selected for model 
validation corresponded to conditions where empirical data 
were available. Model predictions for two additional condi- 
tions were made increasing the engine-to-augmentor inlet 
spacing to one and two engine diameters. No empirical data 
were available for the two latter conditions and, therefore, 
inlet parameters were simulated using empirical data for 
zero engine-to-augmentor entrance spacing. The test condi- 
tions are summarized in Table III. 

Figures 3 through 7 compare predicted axial pressure 
distributions and radial velocity profiles obtained with the 
<p-u and u-v-p computer models. In addition, the available 
empirical data are also plotted on those figures. Pressure 
and velocity were selected for comparison because that was 
the extent of experimental data available. The velocity 
profiles from both computer models were plotted for the grid 
lines closest to the locations of the experimental data. 

In the cases where empirical data were not available, various 
representative velocity profiles were plotted for both 
models. Experimental pressure profiles were available only 
along the top of the augmentor tube. Predicted axial pressure 
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profiles are presented for various radial positions. These 
locations are given as fractions of augmentor tube radius 
(R ) . For example, the pressure distribution labeled 

3 . 

R = 0.96 R indicates that that distribution is along an 
axial grid line located at a distance 96 percent of the aug- 
mentor radius from the axis of symmetry. Two ip-tu model axial 
pressure profiles are depicted for each test condition. One 
profile (R = 0.38 R ) lies in the quiescent flow region 
between the engine and augmentor and was previously found 
by Walters [6] to be the only location which produced a 
reasonable estimate of the measured profiles. The second 
axial pressure distribution (R = 0.13 R ) was located at 
less than one engine radius of the center line. Three u-v-p 
model axial pressure profiles are presented for each test 
condition. One distribution (R = 0.04 R a ) is near the axis 
of symmetry. Another profile (R = 0.28 R ) runs along the 
top of the engine and through the turbulent mixing region. 

The third profile (R = 0.96 R ) is close to the augmentor 
wall. 

2 . Results and Discussion 

a. Test Case I - Low Flow Rate/Zero Spacing 
(1) Velocity Profiles 

As depicted in figure 3, both models pre- 
dicted virtually identical velocity profiles at each station 
along the augmentor tube. There was close agreement between 
the predicted and the experimental velocity profiles at the 
exit of the augmentor. The latter result was expected since 
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the velocity profile has become fully developed near the 
augmentor exit. 

(2) Pressure Profiles 

(a) ip-oj Model - This test condition had 
the lowest flow rates and offered the best chance for agree- 
ment with experiment. The initial pressure drop and the 
pressure rise are underpredicted in the outer flow region 

(R = 0.38 R ) . Nearer the axis of symmetry (R = 0.13 R ) 

a. a. 

the initial pressure decrease is significantly overpredicted, 
but the magnitude (from minimum to maximum) and profile of 
the pressure rise were in good agreement with the experimental 
wall profile. The predicted pressure curves leveled off 
earlier than the experimental data. The difficulty in 
obtaining good pressure profiles with this model are evi- 
dent. Near the augmentor exit there should be negligible 
radial pressure variations. The solution was converged in 
all dependent variables (Table IV) and negligible pressure 
variation was calculated in the radial direction near the 
augmentor exit. However, the pressure profiles along the 
two axial locations presented did not become equal near the 
augmentor exit. The latter resulted from the large errors 
in the predicted profiles near the augmentor inlet. 

(b) u-v-p Model - All three u-v-p pro- 
files predicted an initial pressure drop above the engine 
which was not indicated experimentally. Both the pressure 
drop and rise within the augmentor were underpredicted. 
However, the primitive, variable model more accurately predicted 
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the location of the minimum pressure and seemed to level 
off at approximately the same augmentor position as the 
experimental data. Walters and Netzer [6] have previously 
shown that the model predicts that mixing is nearly complete 
at the location where the pressure profile levels off and 
that, simultaneously, the velocity distribution approaches 
a fully developed profile. Experimental data confirmed this 
characteristic. Figure 3 indicates that, as anticipated, 
the u-v-p model can more accurately predict the location 
at which turbulent mixing is complete. It also predicted 
very little pressure variation with radial augmentor posi- 
tion as is known to be true experimentally. 

b. Test Case IIA - Medium Flow Rate/Zero Spacing 

(1) Velocity Profiles 

The results are presented in figure 4 and, 
again, both computer predictions were very similar and 
agreed with the limited experimental data at the augmentor 
exit plane. 

(2) Pressure Profiles 

(a) Model - The pressure profile 

nearest the center line became more unrealistic for this 
higher flow rate condition. The initial pressure drop was 
greatly exaggerated; however, the magnitude of the pressure 
rise was again in good agreement with experiment. The pres- 
sure profile in the quiescent flow region (R = 0.38 R ) did 

a. 

not agree with the experimental curve. A premature pressure 
drop was predicted and this model underestimated both the 
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pressure drop and rise in the augmentor tube. In addition, 
the minimum pressure point was predicted to occur about one 
engine diameter downstream of the experimental minimum. 

(b) u-v-p Model - As for Case I, the 
u-v-p model predicted a small pressure drop above the engine 
which was not indicated by the experimental measurements. 

It also consistently underpredicted the augmentor pressure 
drop and rise. The slope of the pressure rise, however, was 
in reasonable agreement with experiment. 

c. Test Case IIB - Medium Flow Rate/One Engine 

Diameter Spacing 

(1) Velocity Profiles 

Experimental data were not available for 
this test condition. The predicted profile are presented 
in figure 5 and show that the two computer models predicted 
velocity profiles which were in close agreement. The primi- 
tive variable model, however, predicted a slightly greater 
initial jet spreading at the augmentor entrance and required 
a slightly longer duct length to obtain a fully developed 
profile . 

(2) Pressure Profiles 

(a) ip-w Model - As for Case IIA, the 
centermost pressure distribution greatly exaggerated the 
initial pressure drop. In this case, the pressure rise also 
appeared to be much too rapid in comparison to the R = 0.38 R 

cl 

profile and the u-v-p profiles. The profile in the quiescent 
flow region was in reasonable agreement with the wall profile 
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obtained with the u-v-p model. The slopes of the pressure 
rise and the minimum pressure predicted by both models were 
nearly identical. The \p-u model, however, again predicted 
the minimum pressure to occur somewhat farther downstream 
than did the u-v-p model. 

(b) u-v-p Model - For this case, the 
pressure profile closest to the augmentor wall had a signi- 
ficantly lower minimum pressure than the other profiles. 

d. Test Case IIC - Medium Flow Rate/Two Engine 

Diameter Spacing 

(1) Velocity Profiles 

No empirical data were available for this 
test condition. The computer predictions are presented in 
figure 6. As for Case IIB, the primitive variable model 
predicted greater initial jet spreading. In this case, 
however, the velocity profiles became flat and the 

pressure profiles leveled off considerably upstream of the 
u-v-p predicted profiles. 

(2) Pressure Profiles 

(a) 4i- a) Model - The centermost pressure 
profile was completely unrealistic. The quiescent region 
profile indicated a larger pressure drop and rise and a 
minimum pressure point farther downstream than the other 
model. In addition, the ty-u model profile leveled off much 
earlier. 

(b) u-v-p Model - Again the minimum pressure 
was obtained for the profile closest to the augmentor wall 
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and was located just inside the augmentor entrance. All 
the primitive variable curves indicate a more gentle aug- 
mentor pressure rise, leveling about midway down the augmentor 
tube . 

e. Test Case III - High Flow Rate/Zero Spacing 

(1) Velocity Profile 

As indicated in figure 7, substantially more 
experimental data were available for this test condition. 

In the experiment the nozzle was operated with a pressure 
ratio l ess than critical for one dimensional isen- 

tropic flow. However, using the experimental nozzle flow 
rate and approximating the nozzle flow as one-dimensional and 
isentropic resulted in a nozzle exit Mach number of approxi- 
mately 0.95. This condition was imposed on the models. As 
has been observed for the previous test cases, the velocity 
profiles for both models are quite similar and in reasonably 
good agreement with experiment. However, the predicted pro- 
files for both models tended to flatten a little too rapidly. 
The primitive variable model again predicted slightly more 
initial jet spreading at the augmentor inlet and less mixing 
down the augmentor tube than did the ^-o> model. The experi- 
mental data more nearly agreed with the ip-u model at the 
augmentor inlet and with the u-v-p model downstream. 

(2) Pressure Profiles 

(a) Model - Again the center pressure 

profile was completely unrealistic. The quiescent region 
pressure profile predicted the augmentor pressure drop to 
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begin prior to the engine exit, and underpredicted the 
augmentor pressure drop and rise. The minimum pressure 
position was substantially displaced down the augmentor 
tube and the slope of the pressure rise did not agree with 
experiment . 

(b) u-v-p Model - Aside from the unexplainable 
pressure spike at the augmentor inlet for the centerline 
pressure profile, all three u-v-p pressure profiles were 
quite similar. They underpredicted both the augmentor 
pressure drop and rise. The slope of the pressure rise was 
in good agreement with experimental data within the first 
half of the augmentor tube. Both the \p-u and u-v-p models' 
predicted pressure profiles leveled off before the experi- 
mental curves. 

3 . Comparison of Computational Accuracy and Required 

Computer Time 

The utility of any computer program using numerical 
methods is reflected by the amount of CPU time required and 
the ease of arriving at converged solutions. Table IV 
compares the percentage change in variable magnitude on 
successive iterations and the required CPU time. 

A considerable savings in CPU time was obtained using 
the line-by-line iterative procedure of the primitive varia- 
ble model in lieu of the point-by-point (Gauss-Seidel) method 
employed in the \p-u model. 

At low flow rates, the convergence was quite similar 
for both models. However, at higher flow rates the u-v-p 
model had better convergence in less time. 
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4 . Computer Related Problems 



Both models required that the proper relaxation 
parameters be selected in order to obtain convergence. The 
lack of any procedure for selecting the proper relaxation 
values makes this process quite time consuming. Previous 
research using the \p-ui model facilitated the selection of 
these parameters for that model. It was found that the 
u-v-p model was quite sensitive to the calculated "pressure 
corrections". Obtaining the correct underrelaxation value 
for pressure proved to be the key in arriving at a converged 
solution for the primitive variable model. 

The line-by-line iterative procedure used in the 
u-v-p model was, as expected, quite good in propagating 
disturbances downstream when iterating from left to right. 

A downstream disturbance is propagated upstream by successive 
iterations through the entire field. This fact, at least 
for a geometry incorporating a sudden contraction, made the 
convergence dependent on the number of traverses on each 
radial line. An excessive number of traverses would cause 
divergence. The number of traverses on each line was con- 
trolled in two ways. After each traverse, residual factors 
were calculated for each variable and the largest residual 
factor was compared to a preset value. If the largest 
residual was less than the preset value, the program advanced 
to the next radial line. The program would also advance 
when a preset maximum number of traverses had been completed 
on any radial line. To aid convergence, the program was 
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held to a few traverses on each line until several field 
iterations had allowed the presence of the contraction 
wall to be "felt" upstream. The number of traverses on a 
line was then increased. It was additionally found that 
when working with coaxial flows with radically different 
inlet velocities, the normalizing factors (which were based 
on average inlet conditions and used to calculate the residual 
values on each line) resulted in excessive traverses being 
made in regions of high flow velocity. Repeated calculations 
on radial lines which had already converged often caused 
divergence. Adjusting the normalizing factors downstream 
of the engine exit alleviated this problem. 

The primitive variable model demonstrated some con- 
vergence difficulty in the recirculation region adjacent to 
the sudden contraction. This problem could have resulted 
from the relatively large normalizing factors used in this 
local region of low velocity or, as suggested by Launder 
and Spalding [12] , it could have possibly, been due to the 
inadequacy of the empirical constants (Table I) in the K-e 
model to adequately describe the flow in this quiescent zone. 

As with any finite difference numerical solution, 
grid spacing was found to be critical. To aid convergence, 
the grid spacing approaching any solid boundary was decreased 
and grid lines were packed into regions of large property 
gradients. Care must be exercised in picking successive 
grid sizes. Gosman, et al [1] recommended that for the 
- oj model successive spacing should not increase by more 
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than about a factor of 1.5. This restriction was also 
employed for the primitive variable model. In simulating 
the TJTC, a 30 by 30 grid system was utilized for the primitive 
variable model. A 43 by 40 grid was required for the ip-w 
model . 

5. Summary of Results 

In most cases the predicted velocity profiles for 
both the ip-ui and u-v-p computer models agreed with each 
other and with the available experimental data. The u-v-p 
model predicted more initial jet spreading at the augmentor 
inlet, especially in those cases where the engine exit 
plane was not flush with augmentor inlet. Both models pre- 
dicted that a flat velocity profile was obtained where the 
pressure profile leveled off. For the one case in which 
experimental data were available over the entire length of 
the augmentor tube, the predicted velocity profiles seemed 
to flatten slightly faster than the experimental data. 

In general, the model demonstrated poor pressure 

prediction capability. For low flow rates the pressure rise 
on the wall from minimum to maximum was accurately predicted 
for a pressure profile calculated near the axis of symmetry. 

For higher flow rates the centermost predictions became 
unrealistic although converged solutions were obtained for 
all primary variables. For the predicted pressure profiles 
along axial lines that were in the quiescent flow region, the 
ip-u> model characteristically underpredicted the pressure 
drop and rise along the augmentor wall. The minimum predicted 
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pressure location was typically somewhat downstream of its 
experimentally determined position. In addition, contrary 
to experimental evidence, the model predicted the start 

of the pressure decrease significantly upstream of the engine 
exit plane. 

The primitive variable model showed little radial 
pressure variation with radial augmentor position except 
near the engine exit plane. It consistently predicted an 
erroneous but small pressure drop above the engine. It 
seemed to accurately locate the minimum pressure position 
and to predict the rapid pressure drop at the augmentor 
entrance. It consistently underpredicted the pressure drop 
and rise in the augmentor but tended to accurately predict 
the slopes of the rising pressure profiles. Application of 
the K-e turbulence model with fixed parameters (C , C^, C^/ 
a , o ) to the test cell flow conditions may be the major 
reason for the lack of quantitative accuracy in the predicted 
pressure profiles. Further experimental data (for example, 
turbulence intensity measurements) are needed to determine 
the models applicability to the low velocity region sur- 
rounding the engine. As indicated in figures 3 through 7, 
there were some significant differences between experimental 
and predicted "pressure differentials" at the augmentor exit. 
Calculations indicate, however, a maximum of two percent 
error in absolute pressure at that location. Therefore, 
continuity was satisfied within acceptable limits using the 
virtually identical predicted and experimental velocity profiles. 
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The primitive variable model appears to reasonably 
predict the TJTC flow field up to sonic engine exhaust 
conditions. These predictions include realistic pressure 
distributions and require substantially less computer time 
and fewer grid lines than the ip-ui model. 

B. SOLID FUEL RAMJET 
1 . Introduction 

The purpose of this portion of the study was to 
use a primitive variable, finite-difference computer program 
to determine the flow within a solid fuel ramjet combustor 
with emphasis on the aft mixing chamber. The model was then 
used to investigate the effects of air flow rate through the 
fuel port on the combustion in the aft mixing chamber. As 
previously explained, an aft mixing region allows further 
combustion aft of the fuel grain. This process increases 
combustion efficiency. Lowering the inlet flow rate increases 
the fuel-air ratio within the combustion region. Bypass air 
can then be injected into the aft mixing region. The latter 
procedure can be used to appreciably increase fuel loading. 
Previous work at the Naval Postgraduate School [5,13,14] 
modeled a SFRJ with a computer program utilizing as 

primary variables. Numerical instabilities, however, pre- 
vented the use of the 41 -w model to predict the flow in the 
aft mixing region. The results of that investigation and 
some empirical data were available for comparison with 
the predictions from the primitive variable model. 
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Several factors were anticipated which could 
contribute to differences in the predictions of the two 
models and the empirical data: 

a. Some of the experimental data were measured. in 
cold, non-reacting flows. 

b. The incorporation of the aft mixing chamber into 
the primitive variable model could influence the flow upstream 
in the combustion chamber. 

c. In the model, a wall value of turbulent 

kinetic energy (K) was specified through a slip factor 
such that K w = -1.0 or -0.39*Kp, depending on the magnitude 
of the turbulent Reynold's number. In the u-v-p model the 
boundary condition for K at the near wall node (p) was 
specified in terms of the wall shear stress. In addition, 
in the primitive variable model, the boundary condition for 
stagnation enthalpy at the near wall node was made a function 
of wall shear stress through Reynold's Analogy. These factors 
affect heat flux to the wall and, therefore , the fuel regression 
rate . 

d. The u-v-p model used a 23 by 21 grid in the 
fuel chamber while the ip-w model used a 17 by 25 grid. In 
reality the heat of vaporization of the fuel is a fixed quan- 
tity and, if converged solutions are obtained, the wall heat 
flux should not depend upon the grid spacing. However, it 
has been found [5] that the heat flux to the wall (which is 
calculated using the near-wall grid point) is a function of 
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the grid distance from the wall. This results from the 
assumed behavior of the variables near the wall. The pro- 
cedure employed in this study was to adjust the heat of 
vaporization to match the empirical fuel regression rate 
at one air flow rate and to use that value for all other 
flow rates. If the model is realistic, fuel regression rate 
should then vary with air flow rate in agreement with 
experiment. 

2 . Results and Discussion 
a. Regression Rate 

Figure 8 shows that the fuel regression rate 
predictions of the u-v-p and ij>-u> models are quite similar. 
Both predict the peak regression rate upstream of experiment 
and have similar slopes. This early peak in the regression 
rate results from the model predicting a shorter reattach- 
ment length than was found experimentally [5] . The primitive 
variable model predicted a higher regression rate but, as 
previously discussed, the magnitude can easily be adjusted 
through the value used for the heat of vaporization of the 
fuel . 

Figure 9 shows the effects of increasing inlet 

air mass flux (G = m . / A _ « V. ) on fuel regression rate 

air tp in 

(ff u ) . The regression profile remained the same and, as 
expected, decreased with decreasing G. It has been found 
experimentally [13] that fuel regression rate varies as the 
air mass flux raised to a constant power (r_ G n ) . Boaz 
and Netzer [13] found that this constant was equal to 0.41 
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while Mady, et al [14] found it to be approximately 0.38. 

For the three test cases of this study, the u-v-p model pre- 
dicted the constant, n, to be between 0.31 and 0.34. Thus, 
the primitive variable model appears to correctly predict 
the nature of the change in convective heat flux to the 
fuel surface with air flow rate. 

b. Turbulence Intensity 

Figure 10 compares the predicted centerline 
turbulence intensity (assuming isotropic turbulence) and 
experimental data for non-reacting flow. The primitive 
variable computer model slightly underpredicted the peak 
turbulence intensity while the ip-u> model overpredicted it. 
Both models predicted the peak occurring downstream of 
experiment and both distributions seem to be approaching an 
identical asymptote downstream. The decrease in turbulence 
intensity predicted by the model near the inlet resulted 

from the model over-predicting the velocity increase as the 
air entered the combustor [5] . The u-v-p model overcame 
this difficulty. The differences in the results from the 
two computer models may result from the differences in the 
boundary conditions on turbulent kinetic energy in the com- 
bustor and/or to the effects of the addition of the aft 
mixing chamber on the upstream flow. Both of these phenomena 
have already been discussed. It should also be noted that 
the experimental data used in this comparison was measured 
in a non-reacting flow. 
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Figure 11 shows the effect of decreasing inlet 
air mass flux on turbulence intensity. As anticipated, the 
peak turbulence intensity decreased as inlet axial velocity 
decreased. Each test condition, however, converged on the 
same value downstream. Much additional experimental work 
is required to obtain the turbulence intensities in reacting 
flows; only then can the adequacy of the K-e turbulence 
model be fully evaluted. 

c. Pressure 

Figure 12 shows the effect of inlet velocity on 
the axial pressure distributions for the three primitive 
variable test conditions (Table V) . (The radial scale has 
been expanded to illustrate the pressure variations. The 
maximum pressure variation is approximately 1.2 psi.) Like 
the predicted pressure profiles of the TJTC, the radial 
location of these distributions are given as a fraction of 
the fuel port radius (R^) • As expected, pressure initially 
increased due to jet spreading. This was followed by a 
slight pressure drop as the flow accelerated due to heat 
addition and wall friction. The final pressure rise was 
due to jet spreading in the aft mixing region. 

d. Temperature 

Figure 13 displays radial temperature variations 
in the combustor near the end of the fuel grain and at about 
1.5 aft mixing region diameters down the aft chamber. As 
discussed above, fuel flow rate decreases less than inlet 
air flow rate (r fu ^ G n , n < 1) . Therefore, as air flow rate 
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is decreased, the overall mixture ratio becomes more fuel 
rich, and the developing boundary layer and the fuel layer 
between the diffusion flame and the wall thickens. Thus, 
as shown in figure 13, as the inlet velocity (and, therefore, 
the inlet air mass flux) was decreased, the maximum tempera- 
ture (or "flame") in the combustor moved away from the fuel 
grain and the centerline temperature increased. The maximum 
temperature in the aft mixing chamber was also predicted 
to occur farther from the top wall. 

Figure 14 shows similar data predicted by the 
oj model slightly farther upstream. A significant differ- 
ence between the predictions of the two computer models was 
that the ip-v model predicted a stronger dependence of the 
peak temperature radial location on the inlet air velocity. 

An aft mixing region was not incorporated into the ip-u model. 
Therefore, the boundary layer continued to grow and the point 
of maximum temperature continued to recede from the fuel sur- 
face with increasing axial distance from the initial reattach- 
ment point. The aft mixing region of the u-v-p model caused 
the boundary layer thickness (and, therefore, the location 
of the peak temperature) to become approximately constant 
in the latter portion of the combustion chamber. This was 
the apparent cause of the weaker dependence predicted by 
the u-v-p model of peak temperature location and boundary 
layer thickness on inlet air mass flux. 
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e. Aft Mixing Region 

Figure 15 is an illustration of the predicted 
combustion behavior in the aft mixing region. (The radial 
dimension has been expanded for clarity.) Lines of maximum 
temperature (i.e., the flame sheet location) are presented 
as a function of fuel grain inlet air velocity. It should 
be noted that the aft recirculation zone, which is also 
depicted on this figure, was predicted to be fuel rich and 
did not vary appreciably in size with changing inlet air 
mass flux. As discussed above, the fuel regression rate 
decreased more slowly than inlet air flow rate. Thus, as 
air flux was decreased the mixture entering the aft chamber 
became more fuel rich and the thickness of the fuel layer 
at the end of the fuel grain increased slightly. With high 
air mass flux the flame reached the wall, resulting in 
complete combustion. This condition could be expected to 
produce a high combustion efficiency. For the lower air 
flow rates, the flame did not reach the wall. This would 
result in unburned fuel entering the nozzle and a lower 
combustion efficiency. This effect can also be seen in 
figure 13. As the inlet axial velocity was decreased, the 
maximum temperature region ("flame zone") progressively 
moved farther from the wall. These predictions might be 
used as a first approximation for predicting the "best" 
placement of bypass air dumps in the aft mixing region. To 
predict an optimum location, however, the primitive variable 
model would have to be expanded to three dimensions. 
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3 . Computer Related Problems 



As has been discussed previously, in order to obtain 
results that were in agreement with experiment, the grid 
spacing near the fuel surface was required to be fine and 
the length scale of turbulence was decreased on the combustor 
step face. Because convergence was sensitive to the length- 
to-width ratio of individual control volumes, the small 
radial grid spacing near the fuel surface forced similar 
fine spacing in the axial direction downstream in the aft 
mixing region. A length to width ratio of less than ten to 
one was required. These criteria forced the use of a large 
number of cells, which in turn required a large amount of 
CPU time. A typical primitive variable 40 by 33 grid required 
75 to 80 minutes of CPU time to converge. A typical 
model with a 17 by 25 grid required 35 to 40 minutes of CPU 
time. It must be remembered, however, that numerical 
instabilities prohibited the modeling of an aft mixing 
chamber with the ip-u model. 

The primitive variable model demonstrated some con- 
vergence difficulty in the aft recirculation region. This 
problem seemed to be associated with the continually changing 
velocity profile just prior to the aft expansion (the "inlet" 
conditions for the aft mixing chamber) . This effect was 
suppressed by iterating through the entire flow field several 
times with only a few traverses on each line and then increas- 
ing the number of traverses on the radial grid lines in the 
aft mixing region once the combustor flow field had essentially 
converged. 
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4 . Summary of Results 



In general, the predicted flow fields for the two 
computer models were quite similar within the combustor. As 
anticipated, the primitive variable model allowed the predic- 
tion of the flow within the aft mixing region. This was 
not possible with the ip-u model. The presence of the aft 
mixing region coupled with the few boundary condition differ- 
ences previously mentioned, had some effect on the flow field 
predictions in the ramjet combustion chamber. The most 
noticeable of these were the decrease in dependence of the 
boundary layer thickness and the maximum temperature radial 
location on axial inlet velocity. Many additional empirical 
data are needed to completely assess the validity of the 
primitive variable model in predicting the flow in a SFRJ. 

The primitive variable model, however, reasonably predicted 
the solid fuel ramjet flow field and enabled the simulation 
of an aft mixing region. 
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TABLE II. GOVERNING EQUATION PARAMETERS 
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TABLE III. TUBROJET TEST CELL TEST CONDITIONS 
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TABLE IV. COMPARISIONS OF CPU TIME AND (A<p/(p)-100 FOR LAST 
CALCULATION 



CASE r in (m) m air (kg/ sec) u in (m/sec) G (kg/m -sec) 
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SECONDARY FLOW 
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= .09k / . 03r 

FIGURE 1. TURBOJET TEST CELL GEOMETRY AND BOUNDARY CONDITIONS 
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non-dimensional 

FIGURE 2. SOLID FUEL RAMJET GEOMETRY AND BOUNDARY CONDITIONS 
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FIGURE 3. AUGMENTOR PRESSURE AND VELOCITY PROFILES FOR STRAIGHT 
PIPE INLET AND ZERO ENGINE -AUGMENTOR SPACING 
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FIGURE 4 . AUGMENTOR PRESSURE AND VELOCITY PROFILES FOR STRAIGHT 
PIPE INLET AND ZERO ENGINE-AUGMENTOR SPACING 
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FIGURE 5. AUGMENTOR PRESSURE AND VELOCITY PROFILES FOR STRAIGHT 
PIPE INLET AND ONE DIAMETER ENGINE-AUGMENTOR SPACING 
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FIGURE 7. AUGMENTOR PRESSURE AND VELOCITY PROFILES FOR STRAIGHT 
PIPE INLET AND ZERO ENGINE-AUGMENTOR SPACING 
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FIGURE 8. PLEXIGLASS REGRESSION RATES 
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= 1 78.0 m /sec 
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“ '97.0 m/sec 
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AFT MIXING CHAMBER TOP WALL 
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FIGURE 14. PREDICTED RADIAL TEMPERATURE DISTRIBUTIONS 

(t p-u COMPUTER MODEL, FIGURE 8 OF REFERENCE 5) 
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(PRIMITIVE VARIABLE COMPUTER MODEL) 
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